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ABSIRACT 


Propagation of electromagnetic energy through the atmosphere is a difficult task 
because of temperature fluctuations and index-of-refraction inhomogeneities which de- 
grade the beam's coherence. Understanding this phenomena is of practical importance 
for optical systems. 

This thesis presents an analvtical numerical technique which simulates the effects 
of atmospheric turbulence. The extended Huygens-Fresnel principle was used to simu- 
late wave propagation in a two-dimensional randomly varying medium, which is repres- 
ented by thin, filtered, Gaussian phase screens. The wave optics code implements both 
Fresnel and Fraunhofer propagation, by emploving the fast Fourier transform (FIT) 
meormthm. The analytical spatial coherence length, p, , and normalized intensity vari- 
MINERIN?) 72, of the perturbed electric field. were examined. Results support the concept 
of intensitv saturation for weak scattering cases, however. differences in the values of the 


theoretical and analytical spatial coherence lengths. occurred. 
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I. INTRODUCTION 


Atmospheric turbulence degrades the coherence of electromagnetic energy propa- 
gating through the atmosphere. The refractive-index inhomogeneities associated with 
the turbulent atmosphere induce phase and intensity perturbations across the wavefront. 
Understanding the propagation and scattering of optical waves in random media, 1s es- 
sential for atmospheric laser beam propagation and imaging systems. 

This thesis models the propagation through a random media by means of the ex- 
tended Huvgens-Fresnel principle. A two-dimensional, thin, Gaussian phase screen re- 
presents the randomly varying medium. This wave optics propagation code employs the 
fast Fourier transform (FFT) algorithm in both Fresnel and Fraunhofer forms of prop- 
agation. The Fresnel region incorporates two forms of the diffraction integral, the 
transfer function form for “near field” distances, and the convolution form for “far field” 
regions. 

Previously, numerically intensive simulations of this tvpe required large. super- 
computer systems. However, these computations were performed on a compact, 
desktop computer svstem. 

The model's accuracy was assessed by comparing computer values of the spatial 
coherence length, p, . with values obtained from the analvtical mutual coherence func- 
tion (MCF). Additionally, the concept of intensity variance saturation was examined for 
a single phase screen in the Fraunhofer approximation. The normalized intensity varl- 
ance, c?/[*. approaches a saturation value asymptotically for increasing values of the 
index-of-refraction structure parameter. C2. Theoretical calculations suggest saturation 
for multiple scattering. however, thev sav nothing about a single phase screen realiza- 
tion. This thesis provides results which support saturation effects for single scattering 


Cases. 


II. BACKGROUND 


A. STATISTICAL DESCRIPTION OF TURBULENCE 

1. Random Variables 

Maintaining the coherence of an electromagnetic wave propagaung through the 

atmosphere. requires an understanding of the effects imposed on the wave by the tur- 
bulent medium. The fundamental characteristic of atmospheric turbulence 1s its ran- 
domness, which must be described statistically. In addition to statistical quantities, 
further assumptions must be made of atmospheric turbulence. These include the con- 
cepts of stationarity. homogeneity and isotropy. Stationaritv implies that random 
processes are time independent. Homogeneity assumes invariance under a Galilean 
transformation of coordinates, while isotropic variables are invariant with respect to 
coordinate rotations. [Ref. 1] Mathematically, these two assumptions imply that the 
statistics at two points 7, and 7, depend only on the difference, »,= In — Pali: 

2. Local Homogeneity and Isotropv 

In general. atmospheric random variables do not obey the assumptions of 

stationarity, homogeneity. and isotropy. However. Tatarski [Ref. 2] introduces the 
concept of "local" homogeneous and isotropic random variables. This concept requires 
homogeneity and isotropy within a localized region of size Z4, the outer scale length. 
Furthermore, the difficulties associated. with nonstauonarv random variables; ۸16 ۴ 
moved bv considering random fields with stationary first increments [Ref. 2]. Tatarski's 
method applies to a nonstationary random function whose mean varies slowly with tine, 
bv considering the diflerence of the function at two different locations. The slow func- 
tional chanees do 0۲ 2۱۱۳۲ ۲۳ alle or tiada 

2 


3. Structure Functions 


Tatarski introduces the structure function 
دک‎ 2 P 
Dr, — ту) ھ۶۶٦‎ (1) 


à tensor that is the difference between two quantities. Some important aspects of the 
structure function are; that its general form is valid for any variable, and < > denotes 
an ensemble average taken over all possible point pairs 7, . 7, Assuming homogeneity 
and isotropy. the vector dependence reduces to the magnitude re = ir: ) + ٣ 


structure function becomes 


to 


Dr) = <) -An)X >. (2) 


Kolmogorov showed through dimensional analysis. a simple power law dependence of 


equation (2). over a limited interval called the inertial subrange, as 
G 2 
Dár) = Ä (3) 


Tatarskı introduces the concept of "passive additives”, which are quantities independent 
of position in the vector field, and do not directly influence the dynamics of the turbulent 


medium. Temperature 15 a passive additive and has a structure function of the form 
2 2/3 
Disco: (4) 


As long as r remains Within the inertial subrange, temperature is approximated as a 
passive additive, and equation (4) is valid. Likewise the index-of-refraction is a passive 


cuve With a structure function 


€ 213. E 


Dar) = Cp? (3) 


Exucestheundex of refraction depends on the density of the atmosphere, D, and D, are 


felaved bv 


Dy = Daf HUNE e y (6) 


/ 


4. Covariance, Power Spectral Densities 
In addition to the structure function, other characteristics of random processes 
mande the covanance (or correlation). and Power Spectral Densities. 1t is the interre- 
lationship of these three quantities which provide a useful method for analvzing random 


processes. The covariance between two random variables S and T can be expressed as 
B EIS رر راد‎ (7) 
However, more frequently it is the autocovariance function 
۱۳ ۱۱۳۱ 109)>]>, (5) 


which ıs needed. Furthermore. if T is homogeneous and <T> =Q is assumed, then 


equation ($) simplifies to 


Bi(r) 2 «T(ry)(r)». (9) 


Combining this relation and eguation (2), gives an expression for the relationship be- 


tween the structure function and covariance function as 
Dor) = 2|Втт(0) — Вт). 10) 


In one dimension the covariance function and the power spectral density are transform 


pairs given by 





Too 
W(x) = کي‎ (11) 
and 
do 
Br) = 4L к ОО (12) 


Using the fact that B(r) is an even function and substituting equation (12) into equation 


(10), Dir) becomes, 


„ә 
طے۔ہ‎ 


+00 
Dir) = Ji [1 — cos(xr)] H («)dk. (1 


TOO 


Tatarski develops an expression for the one-dimensional Kolmogorov spectral density, 


elven by 
H(k) 2 0.1224( KUU (14) 
1 


Discussion to this point has been of one-dimensional random processes, how- 
ever these concepts are applicable to three-dimensional cases. Analogous to equation 


(11), Tatarski defines the three-dimensional power spectral density as 


Too Too Too 


Ф(к) = eI BAT, 


+00 +00 +00 


B(F) = ده‎ 





Using the relation 





ST) 


the three-dimensional Kolmogorov power spectral density becomes 
l شک ون‎ / 
۷۶۷ 9 ٦۷ 


B. EM PROPAGATION THROUGH TURBULENCE 
1. Wave Equatior 


(16) 


Based upon Tatarski, Clifford [Ref. 3] develops theoretical results of line-of-sight 


propagation through the atmosphere, directly from Maxwell's equations. Assuming zero 


conductivity, and unit magnetic permeability in the atmosphere, as well as. a sinusoidal 


time dependent electromagnetic field. Maxwell's four equations, in Gaussian units, take 


the form 
m 
Vx<E = ikll. 


— 


VxH - —ikn'E, 


n 


Rewriting equation (22) in the form 
E.Vn! nV. E - Q0, (24) 
and substituting it into equation (23). yields the vector form of the wave equation 
WE + k’n*E + 2V(E -V logn) =0. (25) 


The third term of equation (25) describes the change in polarization of a propagating 
electromagnetic wave. This term 1s negligible as long as the wavelength is small com- 


pared to the refractive inhomogeneities. Thus equation (25) reduces to 
2 一 3 
VE 4 k"n E =0. (26) 


Equation (26) is the vector form of the wave equation describing propagation through 
the turbulent atmosphere. The difficulty in solving this equation lies in the second term 
containing the random variable n. Various methods are available for obtaining solutions 
to equation (26), each of which rehes on several critical approximations. Strohbehn 
[Ref. 4] lists these approximations as: 


1. Negligible depolarization effects. 


{J 


Negligible back-scattering. 
3. Use of the parabolic approximation to the wave equation. 


4. Turbulence ıs uncorrelated in the direction of propagation. 


2. The Method of Small Perturbations-Born Approximation 
Both Tatarski [Ref. 2] and Chíford [Ref. 3]. solve the wave equation in a tussis 
lent atmosphere using the method of small perturbations, which 1s equivalent to the 
Born approximation. This method expands the electric field into a series of decreasing 


amplitudes. and the refractive index into a power series in the form 
E= Egt Eten 23 
n (28) 


Substituting these into equation (26) and equating same order terms, results in two 


equations 


۳ 070 ٦ (28] 





3 E = 0, (30) 


where terms of order mi and higher are ignored. Assuming, as Tatarski does, that the 
unperturbed field is a plane wave propagating in the z-direction represented as 


E, = explikz] , equation (30) becomes 
К БКО = OR کو‎ (31) 


an inhomogeneous partial differential equation with constant coefficients. Its solution 
is the convolution of a plane wave Green’s function with the source term, or inhomo- 


geneous term, given by 


_ exp(idlr — r']) 
dr aa 


PEA 


E(r)- [2^ (F^) exp(ikz!)]. (32) 


ale 
Ja 

5 
3. The Method of Smooth Perturbations-Rytov Approximation 

In addition to the Born approximation, Tatarski develops the Rytov approxi- 


mation Which assumes a solution to equation (30) of the form 
E — exp( Y) 2 exp(X 4 i5). (33) 
or simplv 
Dx اکا‎ (34) 
where A is the amplitude given bv 4 2 exp X. Applying equation (33) to equation (30) 
MINNA TVIdTn£ bv E, vields the Rytov solution grven by. 
V E 4 k^n^(r) s V log E 4- (Vlog EY - K^n'(r). (35) 
Tatarski further shows that both methods of approximation are eguivalent. 


C. HUYGENS-FRESNEL THEORY 

As we have seen, the theories of Tatarski and Clifford use the differential eguation 
approach to solve the problem of propagation through a turbulent medium. However, 
Lutomirski and Yura [Ref. 5], approach this problem in terms of integral equations 
which use an extended Huvgens-Fresnel theory. This technique is equivalent to a dif- 


ferential equation approach, but it is easier to integrate and simulate using FFT tech- 


niques. Lutonurskt and Yura, develop an extended HMuygens-lresnel (icon wees 
introducing a random phase term for turbulence in the Huvgens-Fresnel integral which 
is developed in standard optics texts like [Ref. 6]. This additional phase perturbation 
takes the form of the Rytov approximation, e* . The extended Huvgens-Fresnel integral 


IS, 


—ik | eGklr —x]) 
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27 r—r 





(36) بے میں مد = E(F)‏ 


In the geometrical optics limit, Fermat's Principle is, V fr(z)dz. With a power series 
expansion of e*, equation (36) reduces to the Green's function solution of Tatarski and 
Clifford given in eguation (32). 

Recently, Martin and Flatte [Ref. 7] presented an atmospheric turbulence algorithm 
Which uses the differential equation approach that has as a solution, the extended 
Huvgens-Fresnel Integral. A filtered random Gaussian phase screen introduces the 
phase perturbations while the algorithm’s path integral method, incorporates a multi- 


screen transfer function form of Fresnel propagation. 


D. MUTUAL COHERBNGE'FUN CIO ا‎ 

The effects of atmospheric turbulence can be expressed in terms of two functions, 
the Mutual Coherence Function (MCF) and the Modulation Transter Function ooe 
Lutonurski and Yura [Ref. 5] derive the first concept bv considering the average intensity 
<J(r)> = <E'(F)E0,)> . of eguation (56). What results is an average intensity sS 
a product of the autocovariance of the aperture and the atmospheric NICF. The atmo- 
spheric MCF term has the Rytov form « exp( 1" -- 1)» where the 1" rel Gr SSUG MA 
complex phase factor at the r, coordinate and V ' " corresponds to the»; coordinate mms 
term 1s log-normaliv distributed. as long as. ¥’ 1s composed of Gaussian variables. Using 
this fact and results ın Fried’s work [Ref. $ ], the atmospheric MCF was written in terms 
of the wave structure function D(p ), given by 


RE FELD 
€ 


< (37) 


Eid 
<< 6۸۱ ٠٣۳ 


where D(p)= Dy(p) + D,(p). Lutomirski and Yura apply the structure function for a 


000 qa 


ON 





lie 
D(p) =2.91k PP | C%z)dz, (38) 


to equation (37) which can be written as, 


)39( اوہ 5 (- exp!‏ = ٹہ 


where p, , the lateral cohernece length, given by 


E 
po =|1.464? | Ce)”, (40) 


2 


Mere A =—— aand C? is the index of refraction structure parameter along the optical 
7 

puuneneth- I. represents the distance where the spatial coherence of the wave drops to 

€! point of the MCF. [Ref. $] 


E. MODULATION TRANSFER FUNCTION 

lomiursk)and Yura' s concept of MCF is closely related to the MTF of Fried’s 
Nus The MCF and MTF are in fact the same function but expressed in terms of 
ENrent variables. The MCF is measured in the coordinates of the propagation field 
and has dimensions of distance, while the MTF, is measured in the image plane and has 
the dimensions of spatial frequenev. Both are equivalent under a transformation be- 
tween the two planes by letting p>/Rf where R is the focal length of the optics and f is 
the spatial freguencv. This transformation is valid under the Wiener-Khintchine theo- 
rem since the lens Fourier transforms the incident electric field at the aperture to the 


image plane. Fried's expression for the atmospheric long term MTF is 
AR i 
MTF() = o] EZ FR - (41) 


Eure, — 2.1p,. 
in calculating the MTF, two distinct cases exist. These are a short term and long 
term MTF. The short term exposure describes the evaluation of the wavefront in suffi- 


cientlv short time intervals such that the turbulence appears frozen. The long term 


MTEF, is a single long time integrated exposure, taking into account every turbulence 
configuration. This thesis focuses on the method prescribed by the short term MTF. 
[Ref. S] 

The analvtic distinction between the two cases lies in the manner in which the 
wavefront distortions are handled. Specifically, the distortion attributed to a random tilt 
of the wavefront. Tilt results from the varving phase fluctuations across the aperture 
which accumulate along the optical propagation path. These fluctuations produce im- 
age motion in the focal plane of the receiver. For a very short exposure, the tilt is ex- 
tracted, by fitting a mean square two-dimensional flat plane to the electric field and 
rotating it through an angle so that the mean wavefront is normal to the direction of 
propagation. This introduces a phase shift, resulting in the displacement of each curve 
about the optical axis. Fried's [Ref. 8] development of this theory suggests a higher 
MTF at all spatial frequencies for the short term MTF. The short-exposure MIs 


near-field and far-field cases is respectively given by 


Ana "HR 
«(f)», = toff) exp EX E p 一 i "n 





Fo 
Rf || LIRE 
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where 7, 15 the MTF of a diffraction-limited lens, and the exponential term corresponds 


(42 








to the atmospheric MTF. These equations predict a near field short terny WEP tiem 
start at one, declines to a minimum, then increases to unity, at the optical cut off fre- 


quency. Figure 1 illustrates this phenomena. 


F. DIFFRACTION INTEGRAL 

This section presents the analytical work concerned with the development of my 
propagation algorithms. It begins by illustrating the approximations that were used to 
manipulate the diffraction integral and then proceeds with the analysis required to 
transform the solution of the Helmholtz equation into two forms of the Huygens-Fresnel 
principle. The two forms are, first. a convolution form suited for long distance propa- 
gation and second, as a transfer function form for short distance propagation. Some 
initial assumptions made by Roberts (Ref. 9) include the following: 

I. Light propagates in the k direction. 


2. The wave amplitude is known in the xv-plane at z=0. 
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Figure 1. Short Term Mutual Coherence Function for an 8 x 8 Subaperture. 
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3. Polarization effects are negligible. 


4. The electric field amplitude is a scalar function V(r,z). 


Following Roberts, the analysis begins with the Huygens-Fresnel integral for the 


propagation of light waves given as 





۲, 2( = Jovis. امه‎ El rar | (43) 


where / is the wavelength of light and p is a vector in the aperture plane, r 1s a vector 
in the image plane, and z is the propagation direction. From the Fresnel approximation 
where, r < < zand p < < z, factor z? from the square root and expand by a binomial 


expansion. Equation (43) becomes 
/ 7 2 
Hr.) =" F(p, 0) exp NIU | 1. 
7 Lcd F—D 
lero 0) ae (1 HEY (44) 


ORO exp( E 7-7).‏ اس | esp =E‏ — ا 
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Since the exponential phase factor outside of the integral does not affect intensity 


measurements, equation (44) 1s 


mar E | إمى ره .ماج‎ Lira | (45) 


Expanding the quadratic term gives, 


1 B 1 ja [ave 0) exp| A exp zl gus i | 
AL /DZ jo 7 


This is the convolution form of the Fresnel integral, 3 ي٦‎ 


(46) 








Fraunhofer integral except for the quadratic phase factor in the integral. 
To obtain the transfer function form, the analysis begins with denoting I (f.z) as the 


Fourier ۱۲۵۳۸۹۱۵۲۱۸ ۵ ۱۲2 oi 
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V(f, z) = | Jar exp( —2nif + F)V(F, 2), ۳ 


where F(r, z) is given in equation (45). Interchanging the order of integration yields 








PO, =)= FÉ | арр, |» امه مدمه‎ ۲-2۳ | (48) 
Next, a change of variables is made where 
F'—r—p. (49) 


Substituting this into equation (4$) gives the following equation 





Vif, z) dpV(p. [ur امہ‎ —2zif - (F -- p) ] exp] E E 
B NE (30) 
a (p, 0) exp( — 2zif « jar exp( —2z if « r') exp = ^| 
The r ' integral 1s replaced with its Gaussian transform pair, 
inz exp( —injzf), (31) 
giving 
o D SENP =ni (۷ (p, 0) exp( —nif » 0). (52) 


The inverse Fourier transform of equation (52) yields 


(F, z) = |F expl (2nif ©) exp( —irÃzf "7 (p, 0) exp( —2zif * p). (55) 


This expression is the transfer function form of the diffraction integral, which is equiv- 
alent to the solution of the differential equation approach used by Martin and Flatte 
[Ref. 7]. Reviewing the form of equations (46) and (53), the need for two equations de- 
scribing different propagations. 1s obvious. In one instance, the propagation distance Z 
enters in the denominator of the exponential term. [t is at long distances that the ex- 


ponential term varies slowly. On the other hand. equation (53) is suited for short 


propagation distances, as z enters into the numerator of the 


anons in ths caso 
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IM. NUMERICAL SIMULATION MODEL 


This numerical simulation, modeled wave optics propagation of an electromagnetic 
wave through a random medium, represented by a two- dimensional Gaussian phase 
screens. Techniques in this model required certain “tools” and their testing. These 
"tools" included a need for a reliable random number generator, used in generating the 
random phase screens, and an efficient two-dimensional FFT, used in approximating the 


diffraction integrals. Dut first, a discussion of the experimental arrangement is needed. 


A. EXPERIMENTAL ARRANGEMENT 

Due to the extensive numerical calculations in this simulation, a Compaq desk pro 
80386-20 computer was used. It features a 64 megabyte hard drive and 16 megabytes 
of memory. In addition to the 20-MHz 80387 coprocessor, a Weitek 1167 math 
coprocessor was added to enhance execution speed. The 32 bit Fortran-3$6 compiler 
was from Silicon Vallev Software (SVS) and uses Phar Lap Software to extend the op- 
erating svstem bevond one megabyte. The math and graphics packages were produced 
by Scitech Scientific. The Compay has a 60 x 480 pixel VGA graphics monitor. Both 
Meera? Laser Jet Series 11 and Panasonic KA-P10921 multi-mode printers were used in 


this arrangement. 


B. COMPUTER PRELIMINARIES 
I, Random number generator 

The main purpose of a computer simulation 1s to approximate natura! phe- 
nomena. Jo make things realistic, random number sequences were used to introduce 
stochastic variations. One might ask, what minimal criteria should a particular random 
number generator satisfv? Certainly the most important criterion is that the sequence 
of numbers is sufficiently random. Other criteria are uniformity, reproducibility, mini- 
mum memory, fast. non-repeating, and statistically independent. 

The more difficult characteristic to satisfv is statistical independence. Thus a 
series of tests were needed to provide a quantitative measure of the generator's per- 
formance. There are two kinds of statistical tests: empirical tests and theoretical tests. 
Empirical tests focus on how the computer manipulates groups of numbers from the 


sequence and evaluates certain statistical quantities. Perhaps the best known of al} sta- 


tistical tests are the “Chi-Square” tests. Theoretical tests, on the other hand, establish 
characteristics of a sequence using methods based on recurrence rules. [Ref. 10] 

For the purpose of this simulation. only empirical tests were applied to the 
SVS-Fortran random number generator called Ran/I, . The following five tests were 


considered: 


1. Frequency Test. This test determines whether or not the sequence of numbcr aoii 
uniformly distributed as L(0,1). 


2. Serial Test. This test is an extension of the frequency test tO two dimension» 
matrix form. 


Lagged-Product Test. This test checks for correlations between successive numbers 
Over a Senay period: 


J Rün Pests: 


OI 


a. Runs up and down. This tests for long increasing and decreasing seguences of 
numbers. 


b. Runs above and below the mean. This tests for long seguences with values 

consistently above or below the mean 
The results of the five tests are provided in Appendix A. Of these 11۲6 ۱۳:۳ ۴ 
Lagged-Product and Runs Tests are the most critical when simmulating atmospheric tur- 
bulence. These three tests determine whether or not correlation 15 1111۲ 00۱۱6۶ Is 
random number generator which produces erroneous results in the simulation. The 
SVS-Fortran random number generator met the test criteria and proved to De ORCOS 
better generators. However, this generator has one significant draw back. The random 
sequences begin at one of two different values depending on whether the seed value 1s 
positive or negatiave. Thus it 1s critical that the random numbers be called continuously 
in a loop to avoid restarting the sequence. thereby introducing unwanted correlation. 

2. Fast Fourier Transform 

The most repeatedly used algorithm throughout the numerical simulation Was 
the fast fourier transform. Therefore, 1t was necessarv to use the most efficient algorithm 
available. The Scitech Scientific math package provided several options, with subroutine 
FFT2C . best suited for this numerical simulation. This subroutine uses a complex array 
input. The other FFT considered was a routine coded by Dr. Walters which ۵ ee 
in a demonstration package provided by Infotek. This FFT utilizes real arrays and will 
henceforth, be refered to as suba TM 

Each subroutine was timed for various configurations. Specific timing results 


are contained in Appendix B. Some general results, however. are that subroutine 





FFT2C was faster for one-dimensional cases, with subroutine FFT faster for the two- 
Eimenstonal cases. The decline 1n efficiency of FFT2C can be attributed to the extra 
coding required to Convert between real arrays and complex arrays. Subroutine FT 
was selected over FFT2C since the simulation utilized the two-dimensional form of the 
ГЕТ. 

Other techniques which were emploved to increase the efficiency included in- 
stalling a Weitek coprocessor in the Compag. This reduced the processing time to ap- 
proximately one-third that of the original time. The use of common blocks vice 
dimension statements further decreased processing time by 5%. Finally, a portion of the 
FFT algorithm was modified from wim = sin(ang) , to wim = ./ 1 — wre , with a negli- 
BiBledecline in elliciency by 0.01%. 

Because of the discrete nature of the simulation, there exists problems and lim- 
itations associated with implementating the FFT. One such problem was that of clas- 
sical edge diffraction associated with the phase or amplitude discontinuities at the edges 
of the finite screen. As the propagation distance z increases, the edge diffraction spreads 
toward the center of the screen making more and more of the diffraction pattern erro- 
meous. Buckley |Ref. 11] defines the distance from the ends of the screens where the edge 


diffraction ts important as 
DT (54) 


where z is the propagation distance and Q, 1s the root mean square phase deviation im- 
On the wave by the screen. The severity of edge effects, however, 15 reduced by 
the ahasing introduced in the FFT implementation. Aliasing transforms the linear 
Screen to a “circular” one with the last point associated with the first. This results in a 
continuity of phase and amplitude at the edges of the screen [Ref. 11]. 

The most important limitation. was the finite. spatial range imposed bv the 
maximum available grid size. This places a constraint on the available range of fre- 
КООШО used in the FFT from the lowest given by fa. =1/L . to the highest, 
EL where L is the grid length and n is the number of grid points. Figure 3 on 
page 20 illustrates this setup. The FFT provides a least squares fit of sine and cosine 
fuctions to the perturbed wavefront phases. However, this method prevents an accurate 
representation of the wavefront at low frequencies for a Kkolmogorov x^! ? power spec- 
trum. To alleviate this problem, a subaperture was superimposed at the center of the 


erd. The subaperture helps low frequencies, which contain a large portion of the aimn- 


۳ 


plitude, but hurts the high frequencies bv limiting the inner scale. There exists some high 


frequency edge effects, however, these are minimal. 


C. PROCEDURE 
Figure 2 provides a summary of the coding contained in Appendix C. This con- 
ceptual diagram illustrates an overview of the procedural steps of the Fraunhofer prop- 
agation algorithm. 
l. Input Parameters 
Since this step 1s straightforward, extensive discussion 1s not required. However, 
it is important to note that the input parameters were both fixed and variable. The array 
size and filter value were fixed quantities, but the subaperture size, seed value, C? value 
and propagation distance, took on different values. The actual variable names are doc- 
umented at the beginning of the code. 
2. Aperture Mutual Coherence Function 
The second step in this procedure calculates the aperture MCF. Subroutine 
MCF does this. FIn this subroutine the initial wavefront was represented in the comp m 
asan L x L square array of complex numbers. Centered within this complex eleeuam 
field was a subaperture. The initial complex electric field had a value of zero. every- 
where. except for the real part of the subaperture, which had ۱۳6 ۲۵۱۵ و‎ 
Figure 3 illustrates this. With the electric field created, it was direct Fourier icr ant. oE 
(DIT) by subroutine DFTIFT . The intensiv of the electric field was calculated, and 
then inverse Fourier transformed (IFT). vielding the aperture MCF. 
2 Planai 07 
The complex electric field was created by the same method prescribed in sub- 
routine MCF. It is important to realize that the concept of aperture size Was seem 
two wavs. One way corresponds to the simulated aperture while the other pertains to 
an aperture with physical dimensions. The simulated aperture is actually a matrix or 
grid which directly reflects the dimensioning size. For example, the simular 0٦ 
complex electric field was actually a two-dimensional matrix corresponding to a 256 x 
256 two-dimensional array. From the input parameters, the simulated subaperture takes 
on grid sizes ranging from an $ x $ to 200 x 200 square matrix. 
The second referencing to an aperture refers to an aperture with actual physical 
units. The physical subaperture was assigned a value of 0.3125 meters, which remains 


fixed regardless of the simulated subaperture size. The length L of ٥٣ 
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Figure 2. Conceptual Diagram of the Simulation Coding. 
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field, on the other hand, had variable phvsical lengths depending on the simulated 


supaperture size. The physical value of L was calculated from 





A 
n 
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L-mx—- [meters], ( 
ھ70‎ 
where m is the phvsical subaperture length in meters, mr 1s the integer value of the array 
dimensioning, and isize 1s the integer value for the simulated subaperture. 
4. Phase Screens 
a. Generation 

Phase screens were created in subroutine GGAUS by constructing a L x L 
matrix of Gaussian distributed random numbers. Each position was assigned an inde- 
pendent random number n. thereby, requiring n? random numbers. Since the 
SVS-Fortran random number generator was uniformly distributed, an algorithm pro- 
vided bv Knuth [Ref. 10] transformed the distribution into a Gaussian one. Two inde- 
pendent, two-dimensional real arravs called phaser and phasei were created, which 
represent the real and imaginarv components of a two-dimensional complex phase 
screen. In this algorithm, the imaginary part of the phase screen was set to zero. 

The domain in which the phase screens are created 1s arbitrary, however, 
filtering was done in the Fourier plane. Creating the phase screen in frequency space 
vice real space, reduces the requirement for an additional FFT when transforming from 
real space to frequency space. Martin and Latte [Ref. 7] proposed this method, but it 
creates difficulties in absolute normalization. Further discussion is presented in Chapter 
four. This simulation. on the other hand, generated the random Gaussian phase screen 
۳۰۰۰۰٠٦٠۹٥٠٠. his in turn was DFId to frequencw space where the complex phase 
ET vas filtered and then the filtered phase screen was IF Td. 

b. Filtering 
Phase screens were filtered to obtain the correct power law form. The fil- 


tering function used in subroutine FLTR was, 
Dol) = 22k, (56) 


where 6, 1s the slab thichness and «G5, 2 0.033 Czx-!!5 , which gives the relationship be- 
tween the phase spectrum and refractive-index spectrum [Ref. 2. pp. 101-102.]. Filtering 
was accomplished by multiplving each phase screen spectrum by the square root of 
equation (56). The correct filtering method requires circular frequency filtering instead 


of a linear one, because of two-dimensional isotropy. The correct form is 


K = J? +x, which is radial everywhere except at thé origin, where it 1s zero. ۲ 
This was reflected in the simulation bv setting the position (1,1) equal to zero in each 
two-dimensional real array which makes up the complex phase screen. 

It is important to realize that although the filtering concept is simple and 
straightforward, the actual implementation 1s not. The difficulty arises from the sym- 
metry properties of the digitally filtered phase screen. 

A one-dimensional case offers a simple illustration of this concept. The 
FFT of a complex function which contains only real components, results in 2 6 
function in frequency space about the Nyquist frequency for the real components, and 
an anti-symmetric function for the imaginaries. With these symmetries present in the 
frequency domain, the correct implementation of the filtering is to mimick these sym- 
inetries. Thus a “folding” technique about the Nyquist frequency was required. How- 
ever, when this concept was extended to a two-dimensional case, as in (he Phat sS 
in the simulation, the symmetries imposed by the FFT, were no longer apparent temas 
frequency domain. To simplify the svmmetry requirements, a real phase screen Was 
used. The real and imaginary spectral components were filtered by folding about the 
Nyquist frequency. It was suggested bv both Brigham [Ref. 13] and Martin and Flatte 
that a complex phase screen containing both real and imaginary components. vields two 
entirely distinct phase screens. However, nothing was provided to support this hypoth- 
esis. 

Another important consequence of filtering resides in the units. The phase 
screens were filtered in x units, but the FFT algorithm operates cin و‎ sS 
Therefore, it was necessary to make a change of variables prior to apphving ۵ sS 

he relationship used for thechanse ہہ‎ 275 


The paper by Martin and Flatté [Ref. 7]. specifies an additional normaliza- 

IR 
NA | 
points and A is the sampling interval. Martin and Flatté provided no explanation for 





tion factor of Az! ın equation (20) where A and where N is the number of grid 
the additional Az! in the filtering. It was not included in the filtering code. 
5. Implementation 
After the filtered phase screen was IFT'd into real space, it was introduced into 
the code as a phase screen and multiplied with the electric field. The array phaser. which 
contains the desired phase field, takes the form of the Rytov approximation, e® , in the 
extended Huygens-Fresenel integral. This form assumes that only phase perturbations 


and not amplitude variations, occur. 





6. Propagation Methods 
As indicated in previous sections, the Huvgens-Fresnel technique was used to 
simulate the propagation of light. This was accomplished by applving the FFT to the 
perturbed electric field. Both the “far-field” and “near-field” propagation methods were 
considered. 
a. Fraunhofer 

Of the two different propagation methods, the single screen Fraunhofer 
propagation is bv far the simplest technique to implement and the one implemented in 
ENSNesi. The uniform coherent plane wave at the aperture was FET'd yielding the 
desired diffraction pattern. Looking more closely, one can see that under certain cir- 
cumstances, Fraunhofer propagation 1s just a special case of the long distance convo- 
lution form of Fresnel propagation given by equation (46). There exists two situations 
when this occurs. One is a “far-field” case for large distances, Where the point of obser- 
vation is at infinity. The other case. 1s when a spherical curvature 1s placed on the wave 
at the aperture. This curvature cancels the quadratic phase factor in the Fresnel form, 
at the focal point. 

Both Fresnel and Fraunhofer algorithms were needed for propagation. [t 
was essential to verifv and validate each case before building on the pre-existing codes. 
The Fraunhofer algorithm provided the basis of this simulation. Sinee the Fraunhofer 
diffraction pattern is well-known, it provided a means to verify the existing code by 
comparing the simulated diffraction pattern with theoretical results. Figure 4 and Fig- 
ure 3 illustrate one three-dimensional quadrant of the diffraction pattern of an unper- 
turbed electric field. for two different subaperture sizes. While Figure 6 illustrates one 
three-dimensional quadrant of a perturbed electric field diffraction pattern for a 16 x 16 
subaperture. 

b. Fresnel 

Although the Fresnel propagation forms were implemented but not tested 
in this thesis, discussion is warranted since multi-screen Fresnel propagation 1s predom- 
inantly used in thermal blooming and multiple scattering scenarios. Propagation 
through the turbulent boundary laver also requires Fresnel propagation codes. Two 
different forms are used for Fresnel propagation. which apply a straightforward FFT to 
evaluate them. These two forms however, do not allow for a variable receiving array 
size. In addition, choosing the correet number of sample points 1s essential. An effective 


approach used to resolve these problems is the Fresnel number. 
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Figure 4. Ihree Dimensional Diffraction Pattern for an 8 x 8 Subaperture. 
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Bioure 5. Three Dimensional Diffraction Pattern of a 16 x 16 Subaperture. 
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The Fresnel number is 


TN 
en (57) 
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here Ar is the receiving aperture size, A p is the transmitting aperture size, J is the 
wavelength of light, and z 1s the propagation distance. Implementation of the correct 
Fresnel form depends on the Fresnel number. When the Fresnel number is smaller than 
the number of grid points in the field length, that is, N;<N , the long distance propa- 
Eon algorithm is used. Conversely, when N2N, , the short distance code is 
implemented.[Ref. 14 | 

Ihe long distance propagation code uses the convolution form of the 
diffraction integral. Implementation requires placing a curvature on the electric field 
wavefront at the aperture. The subroutine called quad! , does this. Multiplying the 
phase screen with the electric field. gives a perturbed electric field that propagates the 
Smee cistance by means of one FFI. fn the Fourier plane, the quadratic phase factor 
called quad? scales the electric field. 

The transfer function form, on the other hand, is suited for short distances. 
In this case, the entire propagaton distance 1s divided into equallv spaced slabs. Two 
FFT's are required to propagate the distance of each slab. This is accomplished by the 
following method: 

resi the electric field and phase screen. 


Ex the direct FFT. 


гә 


3. Multiply the field by the propagation transfer function subroutine called ۲ ۰ 


EN Dpiv the inverse FFT. 


This procedure is repeated for each phase screen until the observing plane is reached. 
۱۶۰۲٢ 
Both forms of Fresnel propagation were included in the simulation code 
contained in Appendix C, however, neither form of Fresnel propagation was exercised. 
7. Atmospheric Mutual Coherence Function 

The final step in this simulation calculated the atmospheric MCF. The same 
procedure used to calculate the aperture MCF, was applied to the perturbed electric 
field. with one exception. The difference is that division of the composite atmospheric- 
N MCT bx the aperture MCF, gives the atmospheric MCF. 


IV: "RESULTS 


The main objectives of this thesis were to demonstrate that the simulation gives re- 
sults that provide insight to weak scattering wave propagation and examine the accuracy 
of the simulation code. The limitations imposed by the computer mechanics and the 
actual method of implementing turbulence theory are discussed. Specific areas of inter- 


est include the filter function, MCF, and saturation effects. 


A. FILTERING 

As previously mentioned in chapter three, Martin and Flatté proposed a very dif- 
ferent approach to creating the random Gaussian phase screen. Their method suggests 
creating the phase screen in frequency space, vice real space, to reduce the TEQUIAAN 
number of FFT's from two to one. Since the domain in which the phase screen 1s gen- 
erated, is arbitrary, this approach seems plausible. However. this method proved to be 
awkward. As the subaperture size was successively doubled, it was necessary to increase 
the strength of turbulence, C?, bv a factor of ten each time, in order to producemmm 
identical level of turbulence as in the previous phase screen. Additionally, since the 
phase screen was created in frequency space. and no symmietries were present, ۵ 
tering technique did not reflect any folding about the Nyquist frequency. [tis not eltää 
that ahasing was accounted for in the Martin and Platte algorithm. Hence, the IFT of 
the phase screen appears to result in a statistically incorrect phase screen. Additionally, 
itis not clear how, with one FFT. Martin and Flatté handled ٥٦ 
ization requirements. With a round trip of FFT's, the normalization problems are au- 
tomatically handled. The weakly filtered phase screen, in turn, led to MCF curves which 
erossly overestimated the coherence length. p,. These problems obtained from Martin 
and Flatté ' s approach to simulating turbulence led to the current coding which created 
one phase screen in real space and applied a folding about the Nyquist frequener, in the 
filtering. to account for the svmmetries introduced by the FFT's. 

Kolmogorov theory of atmospheric turbulence predicts that the graph of the phase 
screen power spectral density versus x yield a slope of —11/3. Figure 7 shows that 
equation (356) produced a filtered phase screen that reflects the —11/5 slope, as well as, 
the correct folding technique. Identical slopes were expected for all subapertures, as well 
as, all possible angles which reflect the circular filtering. Other subaperture profiles show 


a consistent slope value of —11/3. Figure $ corresponds to a 52 x 32 subaperture ata 


to 
AI 





45 degree angle while Figure 9 is representative of a 64 x 64 subaperture at a 90 degree 
angle. Isotropy is apparent in that the circular filtering was implemented correctly 
within the error introduced by the randominess in the screen and the ability to linearly 


fit a line through the data points. 


B. MCF 

The MCF of the electric field provides one method of analvzing the accuracy of the 
simulation. To verify that the WCF was correctly computed, the aperture MCF of an 
unperturbed electric field was calculated. This was easily accomplished since the image 
intensity and aperture MCF are transform pairs and are analvtucal for simple square and 
circular apertures. The aperture MCF 1s just the autocorrelation of the aperture func- 
tion Which is evaluated by calculating the area of overlap of two identical apertures as 
they are moved laterally apart. For a square subaperture, the autocorrelation vields the 
triangle function, with maximum value of 1.0 and minimum value of 0.0 corresponding 
est to the subaperture sıze. Figure 10 illustrates the MCF, or autocorrelation. of 
an $ x $ and 16 x 16 square subaperture. 

The MCF corresponding to the atmospheric turbulence was determined bv dividing 
few © of the perturbed electric field bx the aperture MCF. The value of p, the spa- 
Mal coherence length, was determined froni the turbulence MCF curve. Figure 11 11- 
lustrates a simulated p, value of 3.20 mm for a 64 x 64 subaperture with G = Ix10-5. 
The theoretical value calculated from equation (40) vields p, 2 2.31oun.. Although the 
64 x 64 subaperture gave accurate results other subaperture configurations did not. 
EXC. increased, the MCF curves fell off rapidly towards zero for all subapertures. 
Figure 12 illustrates this for a 64 x 64 subaperture. All subaperture configurations were 
run for two differenct C? values. The simulated p, values were plotted against the cor- 
I onumnssubapertures for each case. Figure 13 reflects the C2 2 1x10-%p, values, and 
eee 1-1 is for C2 = 1x10-'*. The desired trend is for the simulated p, values to approach 
Nheoreuca! value as the subaperture size increased. This trend is visible in 
Figure 13 where C? 2 1x10-? and p, = 2.4nun. However, Figure 14 shows continuously 
decreasing p, Values past the theoretical value of 9.6mm. The results of figures thirteen 
and fourteen point to a problem that may involve edge effects or aliasing. resulting from 
undersampling. The p, values. which were on the order of millimeters, were smaller than 
the tens of centimeter distances of the subaperture mesh size. 

In an attempt to pinpoint the problem. the simulation code was changed to increase 


m xuumber of frequencies in the subaperture bv usine a 512 x 512 array. The results 
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Filtered Phase Screen of a 16 x 16 Subaperture. 
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Figure 8. Filtered Phase Screen of a 32 x 32 Subaperture. 
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Figure 14. Subaperture Coherence Lengths for C? = 1x10-". 


were identical to those for the 256 x 256 dimensioning, within the arithmetic error of the 
algorithm. The information provided by the 512 dimensioning was that the inaccuracy 
of the code was not due to an inner scale problem, however. a low frequency, outer scale 


problem may still exist. 


C. SATURATION 

The last area of investigation was whether or not the simulation predicts the satu- 
ration of intensity for an extended medium modeled by a single phase screen realization. 
Saturation is generally considered to be caused by multiple scattering. a scattered wave 
interfering with a distorted wave. One would expect from the Rytov approximation of 
turbulence theory, that saturation will occur even in Fraunhofer propagation. This as- 
sumption stems from the representation of turbulence in the form. e?, which has a 
magnitude bounded between plus and minus one. Figure 15 illustrates that the nor- 
malized intensity variance saturates With increasing turbulence. Theoreticallv. a nor- 
malized variance of one is expected for Rayleigh statistics. However, it is premature fo 
assume that saturation is inherent in Fraunhofer cases, until the MCF results Quem 
fied. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


This thesis simulated the propagation of plane waves through an extended two- 
dimensional random media using a single phase screen technique. The atmospheric 
turbulence had a Kolmogorov x73 pure power law, while the propagation was strictly 
Fraunhofer. Limitations of its applicability were, primarily, the finite spatial range im- 
posed by the available grid size. Diffraction patterns, correlation functions and intensity 
variance saturation at the observation plane, were investigated. 

The results provided bv the simulation suggest general agreement with the turbu- 
lence theorv. Saturation for weak scattering was supported by this model. The MCF 
curves, although not completely correct, provided insight to the theory and illustrated 
problems which are still present in the current coding. Some specific problems include, 
potential errors in the implementation of the filtering technique from edge effects, alias- 
ing and inner scale problems. or from incorrect normalization. 

Anv further research on this topic should begin with resolving the inaccuracies still 
present in the simulation coding. Several possible reasons were presented, however, an 
error in the filtering of the phase screen seems to be the most likely cause bun 
testing can be conducted on the phase screen. to include calculating the phase screen 
variance, as well as. its structure function D,. Bv comparing the simulated phase screen 
with the theoretical structure function for D,. this technique Will verify whether or not 
the simulated phase screen accurately represents turbulence. 

An assumption was made, that aliasing was not a problem. since the problems as- 
sociated with undersampling were not apparent. Aliasing can be tested bv using finer 
grid sizes and observing changes induced by the higher spatial frequencies. 

After the coding 1s Working correctly, both the convolution and transfer function 
form of Fresnel propagation can be implemented and exercised. In addition, the array 
sizes should be modified to incorporate a dimension of 512 or 1024. Finally. during the 
phase screen generation, phasei should be filled with random numbers to give two usable 
phase screens. It is not obvious whether two independent phase screens will be 
produced, or whether the total energy will be distributed among the two phase screens. 
Therefore, testing should be conducted to ensure that each usable phase sereen possesses 


the correct statistics. 
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APPENDIN A. RANDOM NUMBER GENERATOR TEST DATA 


The following table provides the statistical results of the five empirical tests run on 
the random number generator, Ran/I;. The N? values were determined from the Chi- 
Sguare table in Bevington [Ref 15 J. 

The results of the Lagged Product test correspond respectively to the theoretical 
Mean. u, . Calculated mean, 4 , theoretical standard deviation, o; , and calculated 


standard deviation, o. A lag of three was tested. 


Table 1. RANDOM NUMBER GENERATOR TEST DATA 


um po سس ا‎ 


Above and Below the 
Man 


TEST 
Lagged Product 0.250 0.100 0.092 
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APPENDIX B. FFT TIME SERIES DATA 


This appendix contains the results of the comparison between the two FFT sub- 
routines, FF 72C and subroutine FFT. A function routine called SCN DS, was imple- 


mented in the code to provide accurate time measurements. 


Table 2. FFT TIME SERIES DATA 
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APPENDIX C. SIMULATION CODE 


The simulation code contained in this appendix incorporates Fraunhofer and both 


forms of Fresnel propagation. This thesis only exercised the Fraunhofer propagation. 
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THIS CODE PROVIDES A QUALITATIVE VIEW OF BOTH FRAUNHOFER AND 
FRESNEL DIFFRACTION BY OBSERVING THE PERTURBATION IMPOSED ON A 
MONOCHROMATIC PLANE WAVE PROPAGATING THROUGH A TURDULENT 
MEDIUM. THE TURBULENT MEDIUM IS INTRODUCED IN THE FORM OF A 
STOCHASTIC PHASE SCREEN. PROPAGATION OF THE ELECTRIC FIELD IS 
ACCOMPLISHED THROUGH FFT'S. 
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GLOSSARY OF VARIABLE NAMES: 
1. ARRAYS: 


RE - ONE DIMENSION REAL ARRAY OF DIMENSION NR, WHICH IS USED TO 
MANIPULATE THE REAL PART OF THE COMPLEX ELECTRIC FIELD IN 
THE FFT SUBROUTINE. THIS ARRAY TS REPEATEDLY USED THROUGHOUT 
THE CODE. 


RIM- ONE DIMENSION REAL ARRAY OF DIMENSION NR, WHICH IS USED 
TOSMANIPOLATE THE IMAGINARY PART OF THE COMPLEX ELECTRIC 
ШОО MIEN T SUBROUTINE. -THIS ARRAY IS REPEATEDLY USED 
THROUGHOUT THE CODE. 


FIELDR - TWO DIMENSION REAL ARRAY OF DIMENSION NR X NR 
CONTAINING THE REAL PART OF THE COMPLEX ELECTRIC 
FIELD. THIS ARRAY IS REPEATEDLY USED THROUGHOUT THE CODE 


FIELDI - TWO DIMENSION REAL ARRAY OF DIMENSION NR X NR 
CONTAINING THE IMAGINARY PART OF THE COMPLEX ELECTRIC 
FIELD. THIS ARRAY IS REPEATEDLY USED THROUGHOUT THE CODE 


FILL - TWO DIMENSION REAL ARRAY OF DIMENSION NR X NR USED AS A 
DUMMY ARRAY IN THE IFT OF THE POWER SPECTRUM WHICH YIELDS 
ШШЕ ۶۰ 


FIELDM - TWO DIMENSION ARRAY OF DIMENSION NR X NR REPRESENTING 
MEI TAS IIUPENOP THE PERTURBED ELECTRIC FIELD. 


FMCF - TwO DIMENSION ARRAY OF DIMENSION NR X NR REPRESENTING THE 
Peo ee eer ERIURBED ELECTRIC FIELD. THIS ARRAY IS 
ESED INSPEJERMINING. THE. MEF. 


FNORM - TWO DIMENSION ARRAY OF DIMENSION NR X NR REPRESENTING THE 
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INTENSITY OF THE UNPERTURBED ELECTRIC ۶ ٣,۹ 
IS ALSO USED IN DETERMINING THE MCF. 


PHASER - TWO DIMENSION REAL ARRAY OF DIMENSION NR X NR CONTAINING 
THE REAL PART OF THE RANDOM COMPLEX PHASE SCREEN. 


PHASEI - TWO DIMENSION REAL ARRAY OF DIMENSION NR X NR CONTAINING 
THE IMAGINARY PART OF THE RANDOM COMPLEX PHASE SCREEN. 


FMAG - ONE DIMENSION SLICE OF THE PERTURBED MCF USED IN THE 
GRAPHICS ROUTINE. REAL ARRAY. 


DIST = ONE DIMENSION REAL ARRAY REPRESENTING ” 57 
CORRESPONDING TO THE VALUES IN THE ARRAY FMAG. 


VARAIBLES: 


NR - DIMENSION OF THE ARRAYS EXACTLY AS SPECIFIED IN THE DIMENSION 
STATEMENTS IN THE CALLING PROGRAM. INPUT INTEGER. INDICE. 


NZ - ONE HALF OF NES INTEGER 
M - POWER OF 2 IN THE DIMENSIONING. USED IN THE ITI SU ۶ ٦ 


ISIZE - INTEGER VALUE CORRESPONDING TO A PARTICULAR CHOICE FOR AN 
APERTURE SIZE.INPUT VARIABLE. 


NSIZE - INTEGER VALUE CORRESPONDING TO THE ACTUAL APERTURE STZES 


DELMSH = REAL VALUE REPRESENTING THE فان ند‎ ۱ 
PARTICULAR APERTURE SIZE. 


SEED - REAL INPUT VARIABLE USED TO BEGIN A RANDOM SEOBBNCE OB 
NUMBERS FOR THE SUBROUTINE GGAUS. 


NYES - INPUT INTEGER VARIABLE WHICH SELECTS WHETHER TCORBULEN CEs 
INTRODUCED IN THE CODE. 


FILTER = FIXED INPUT VALUE USED IN THE FILTERING FUNCTION: RES 

CN2 - REAL INPUT VARIABLE REPRESENTING THE INDEX-OF -KEPRAG Raa 
STRUCTURE PARAMETER. DETERMINES THE AMOUNT OF TURBULENCE 
INTRODUCED IN THE FILTERING FUNC Tie ۵٥ 

DREC - REAL INPUT VARIABLE REPRESENTING THE نابات‎ ۷ ۱۱۵ ۳ ۱۳۱ 

DTRNS - FIAED VALUE FOR THE TRANSMITTING ۱ IT :+ 1-۵ 


4 - REAL INPUT VARIABLE REPRESENTING THE TOTAL PROPAGATION 
DISTANCE OF THE ELECTRICAS TE ۳۷ 


DELX - REAL VARIABLE FOR THE PROPAGATION DISTANCE TO EACH SLAB. 
NUMSCR - INTEGER INPUT VARIABLE USED IN ۱۳۰۰۲ ۱۳ ا پہ‎ 


PROPAGATION. THIS VARIABLE CORRESPONDS THE THE NUMBER 
Or EQUALLY SPACED SLABS WHICH MAKE UP THE TOTAL DISTANCE. 
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WVL - FIXED VALUE FOR THE WAVELENGTH. 


[ ان‎ < INTEGER VALUE WHICH DETERMINES WHETHER OR NOT THE FIELD HAS 
PROPAGATED THE TOTAL DISTANCE 2. USED IN NEAR-FIELD 
FRESNEL PROPAGATION. 


DX ELDER OF PIS 
0۱۶۰۰" ۷ 1651777۰ 
MODE - REAL VALUE WHICH DETERMINES THE FORM OF PROPAGATION. 


SIGN - REAL VALUE EITHER 1.0 OR -1.0 WHICH DETERMINES WHETHER THE 
۰۰ط‎ ۰۷۰٠ ۷٣٠۰ ٠٠٠۱٢۰٠۱٥۷۷۷٣۹۴۶ ۶۰ IT ALSOFDETERMINES WHETHER 
NORMALIZATION OCCURS. 


FLDM - MAXINUM VALUE IN THE PERTURBED MCF ARRAY. 
FMAX - MAXIMUM VALUE IN THE PERTURBED ELECTRIC FIELD ARRAY. 
GRAPHICS: 


THE FOLLOWING VARIABLE NAMES ARE EITHER SPECIFIC TO THE SVS 
GRAPHICS ROUTINE OR USED TO MANIPULATE DATA FOR GRAPHING: 


MM TRES ANS: GETC, IUNITP, IUNITV, ISYMB, ITNO, MON, NPRIN, 
MISIN ENIN ONE, TITLE, IX, IY, XMAX, ICOLOR, INCR, NTOT. 
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COMMON /BLK1/ RE(256),RIM(256) 

COMMON /BLK2/ FIELDR(256,256),FIELDI(256,256) 
COMMON /BLK3/ PHASER(256,256),PHASEI(256,256) 
DIMENSION NDEX(20) ,FIELDM( 256,256) ,FMCF( 256,256) , FNORM( 256,256) 
DIMENSION FMAG( 130) ,DIST( 130) ,FILL( 256,256) 
INTEGER*2 IREG(9) 

DOUBLE PRECISION PI 

MN PSS NIST 85845555 ,4,4,2,2,3,3,1,1,0,0/ 
DATA RE/256*:0. 0/ ,RIM/256*0. 0/ 

DATA FMAG/130*:0. 0/ ,DIST/130*0. 0/ 

CHARACTER*4 ONE 

CHARACTER*21 TITLE 

CHARACTER*2 ANS,GETC 

DAMA TITLE/ THE SEED VALUE Is= '/ 

DATA ONE/' '/ 

DATA IUNITP/10/ , IUNITV/20/ 
OPEN(4,FILE='LN.DAT' ,STATUS="!NEW' ) 

ISYMB=22 

ITNO=5 

MON=18 

PI=3. 141592653589792 

NPRIN=0 

MODE=0 


DETENER 


999 CONTINUE 


THIS SECTION OF THE PROGRAM SETS UP THE INPUT PARAMETRO 
SIMULATION. 


ve 7 ! 
*)' HELLO... LET US BEGIN THIS SIMULATION BY ENTERING ' 
WRITE(*,*)' SEVERAL INPUT PARAMETER VALUES. 


WRITEC*, 
WRITE(*,*)' 


WRITE(* 


WRITE" W THE VARIABLE WHICH DIMENSIONS THE ARRAY SIZE IS ۱۷ 


WRITE(*,*)'FIRST VALUE TO ENTER. INPUT THE INTEGER VALUE. 
READ * )NR 

WRITE(%*,%)' 
WRITE(* ee THE SECOND VARIABLE OF INTEREST IS "NSIZE".  THIS' 
WRITE(*,*)'VARIABLE DIMENSIONS THE SIZE OF THE PLANAR ELECTRIC' 
WRITE(* PD FIELD. SELECT ONE OF THE FOLLOWING. ' 

WRITE(*,*)' 1. ۳۷پ‎ 

WRITES 2. FOR 64 X 64 

WRITE(*,*)' 3. FOR 32X2" 

7 ٤٥ 4. FOR 16 X 16' 

WRITE(* 20 5. ۳۳ 7 

WRITE SD 6. FOR 4 X 4' 

WRITE Ty 7. FOR ا‎ 

WRITECH,: zs 


READ(**,"*)ISIZE 
IF(ISIZE.EQ.1) THEN 
NSIZE=50 

DELMSH=. 0031 

INCR=3 
ELSEIF(ISIZE.EQ.2) THEN 
NONEZE- 3 

DELMSH=. 0049 

INCRE! 
ELSEIF(ISIZE. EQ. 3) THEN 
NSIZE=16 

DELMSH=. 0098 

HNGR=2 
ELSEIF(ISIZE.EQ.4) THEN 
NSIZE=8 

DELMSH=. 019531 

INCR=2 
ELSEIF(ISIZE.EQ.5) THEN 
NSIZESA 

DELMSH=. 039063 

TNGR=i 
ELSEIF(ISIZE.EQ.6) THEN 
NSIZE-2 

DELMSH=. 0781 

RR 
ELSE 

NSIZE=1 

DELMSH=. 1563 

Pep 
ENDIF 
WRITE i 
WRITE(*,*)'ANOTHER INPUT PARAMETER IS THE SEED VALUE OF THE’ 
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OO OOO OOO O 


رہ دہ رہ 


WRITE(*,*)' RANDOM NUMBER GENERATOR.  INPUT THE SEED VALUE OF +1.0' 
WRITE(*,*)'OR -1.0' 

READ(*,*)SEED 

WRITE(*,*)' 

WRITE(*,*)'FINALLY INPUT AN INTEGER OF VALUE 1 FOR TURBULENCE, OR' 
WRITE(*,*)'0 FOR NO TURBULENCE. ' 

READ(*,*)NYES 

IF(NYES. EQ. 1) THEN 

WRITE(%*,%*)' 

WRITE(*,*)' INPUT THE VALUE OF FILTER' 

READ(*,*)FILTER 

ENDIF 

WRITE(*,%*)' | 

WRITE(*,*)' INPUT THE VALUE FOR CN2' 

READ(* ,*)CN2 
WRITE(*,*)' ۱ 
WRITE(*,*)' INPUT THE RECEIVING FIELD SIZE IN METERS' 
READ(*,*)DREC 
WRITE(*,%)' i 
WRITE(*,*)' INPUT THE TOTAL PROPAGATION DISTANCE IN METERS' 
READ(%* (7 
WRITE(*,*)' 
WRITE(*,*)' INPUT THE NUMBER OF SCREENS TO BE USED EITHER FOR' 
WRITE(*,*)'FRESNEL OR FRAUNHOFER ' 

READ(* ,**)NUMSCR 

IF(SEED. GE. 1.0) THEN 

ONE( 1: 2)='41' 

ELSE 

Wee 1) 2)= =1' 

ENDIF 

1191 20: 21)=ONE( 1: 2) 

N2=NR/2 
WVL=. 5E-6 
M-ALOG(REAL(NR))/ALOG(2. ) 

DERNS=. 3125 

0۲٦۰0 


THE FOLLOWING STATEMENT DETERMINES WHICH FORM OF FRESNEL IS TO BE 
USED 


MODE-INT((2*DTRNS*DREC)/(WVL*2)) 
THE FOLLOWING SUBROUTINE CALLED MCF DETERMINES THE AUTOCORRELATION 
ORSETEEZAPERTURE WHICH IS TO BE USED IN DETERMINING THE ATMOSPHERIC 
COMEKENCE LENGTH. 
CALL MCF(FNORM,NR,M,PI,NSIZE,N2,DELMSH) 
THIS SECTION CREATES THE PLANAR ELECTRIC FIELD 
DO 40 I=1,NR 
DO 40 J=1,NR 


PUELOR I, J)=0. 0 
PIECDI(I,J)=0. 0 


40 CONTINUE 


POPs i I=N2-NSIZE+1] ,N2+NSIZE 
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DO 4l J=N2-NSIZEFI 07 
EIELDR ٦ 
41 CONTINUE 


THIS SECTION DETERMINES THE SLAB THICKNESS(ES) FOR WHICH THE 
ELECTRIC FIELD IS PROPAGATED THROUGH. THIS TS VABTIDZEURZERTE 
FORMS OF PROPAGATION. 


IFCNR. LE. HODE THEN 
DELX=Z/NUMSCR 

ELSE 

DELX=Z 


THE FOLLOWING SUBROUTINE PLACES A CURVATURE ON THE WAVEFRONT TO BE 
USED FOR THE CONVOLUTION FORM OF FRESNEL PROPAGATION. 


CALL QUADI(NR,PI,DELMSH,DX,DY,WVL,DELX) 
ENDIF 


۳ زی رق رق كت كت‎ 62 6٥ن‎ 6(7 C2 @ ۳ GEO 


1000 CONTINUE 


THIS SECTION CALLS OUT GAUSSIAN RANDOM NUMBERS FOR THE REAL AND 
IMAGINARY PHASE ARRAYS. 


DEDM 


CALL GGAUS(NR,SEED) 


THIS BEGINS THE FILTERING PROCESS OF THE PHASE SCREEN 


OS® 


CALL FLTR(NR,DELMSH,CN2,DELX,WVL,FILTER) 


HERE THE INDIRECT TRANSFORM IS BEING APPLIED TO TIEZEPETERE 
PHASE SCREENS. 


رق ری ۳ 


CALL IFTSCR(NR,M,DELMSH) 


THE FOLLOWING ARRAYS USED IN THE FFT ROUTINES ARE ZEROED OUT TO 
ENSURE THAT UNWANTED VALUES ARE NOT LEFT IN THE ARRAYS. 


ODDS 


DO 990 I=1,NR 

RE(1)=0.0 

RIM(T)=0.0 
990 CONTINUE 


THIS SECTION DOES THE ALGEBRA NEEDED TO MESH THE PHASE SCREEN 
TOGETHER WITH THE ELEĠIRIOSETEJED 


exc c) 


DO 50 I=1,NR 
DO 50 J=1,NR 
XA-COS(PHASER(I,J))*FIELDR(I,J) 
XBECOS( PHASER( I , J) )*FIELDI(I, J) 
XCSSIN(PHASER(I, J) )*FIELDR(I,J) 
XD=S IN( PHASER(1,J))*FIELDI(I,J) 
FIELDR(I,J)sXA-XD 
FIELDI(I,J)=XB+XC 
50 CONTINUE 
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(oc cy c) 


PARO ASS Mo oi iO ME DO O 


6 ۵ 


CO 3 


141 


888 


910 


HERE THE FAST FOURIER TRANSFORM IS BEING APPLIED TO THE PERTURBED 
ELECTRIC FIELD. FOR A DIRECT TRANSFORM SIGN=-1.0, AND INDIRECT 
TRANSFORM SIGN=+1. 0. 


SIGN=-1. 0 
CALL DFTIFT(NR,M,SIGN,DELMSH) 


DO 141 I=1,NR 
DO 141 J=1,NR 
REEDIN) -PTEEDR(OT,J)/(NR*DELMSH) 
FIELDI(I,J)-FIELDI(CI,J)/(NR*DELMSH) 
CONTINUE 


THIS PORTION OF THE IF STATEMENT CORRESPONDS TO THE IMPLEMENTATION 
OE تا‎ TRANSFER FUNCTION FORM OF FRESNEL PROPAGATION. THE 
SUBROUTINE CALLED TRNSFER APPLIES A QUADRATIC TO THE FIELD. 


IF(NR. LE. MODE) THEN 

CALL TRNSFR(NR,PI,DX,DY,WVL,DELX,DTRNS) 
ICNT=ICNT+1 

SIGN=+1. 0 

CALL DFTIFT(NR,M,SIGN,DELMSH) 

GO TO 888 

ELSE 


THE SUBROUTINE CALLED QUAD2 PUTS THE DIFFRACTION PATTERN IN REAL 
SPACE COORDINATES. 


CALL QUAD2(NR,PI,DX,DY,WVL,DELX) 
ENDIF 


CONTINUE 


۱۱۱ ۰ OG LOOP DETERMINES THE POWER SPECTRAL DENSITY AND SETS IT UP 
BURVAN ETT TO DETERMINE THE MCF. 


DO 152 I=1,NR 

DO 152 J=1,NR 
FMCF(1,J)=FIELDR(1,J)**2+FIELDI(I,J)**2 
EE J) =0. O 

CONTINUE 


ONCE AGAIN THE ARRAYS ARE CLEARED OF STRAY VALUES 

DO 910 I=1,NR 

RE(1)=0.0 

RIMC1)=0. 0 

CONTINUE 

PIN VERSE EET IS APPLIED TO THE POWER SPECTRAL DENSITY 
SIGN=+1.0 


DO 901 I=1,NR 
DO 911 J=1,NR 
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911 


91 
901 


94] 


951 
793 


69 


29 
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RE(J)=FMCF(1,J) 
RIM(J)=FILL(1,J) 
CONTINUE 
CALL FFT(M,SIGN,DELMSH) 
DO 921 J21,NR 
FMCF(1,J)SRE(J) 
FILL(1,J)=RIM(J) 
CONTINUE 
CONTINUE 
DO 931 J=1,NR 
DO 941 I=1,NR 
RE(I)=FMCF(1,J) 
RIM(I)-FILL(1,J) 
CONTINUE 
CALL FFT(M,SIGN,DELMSH) 
DO 951 I=1,NR 
FMCF(1,J)=RE(1) 
FILL(I,J)=RIM(I) 
CONTINUE 
CONTINUE 


THIS SECTION DETERMINES THE MAXIMUM VALUE AND NORMALIZES THE MCF 


FLDM=0. 0 

DO 89 I=1,NR 

DO 89 J=1,NR 
XMG=FMCF(I,J) 
IF(XMG. GT. FLDM) THEN 
FLDM-XMG 

ENDIF 

CONTINUE 

WRITE(*,*) FLDM 

PAUSE 


THE MCF IS NORMALIZED SO THE MAX: VALUE Sissi ۵ 


DO 29 I=1,N2 

DO 29 J=1,N2 
FMCF(1,J)=FMCF(1,J)/FLDM 

CONTINUE 


THE ATMOSPHERIC MCF IS DETERMINED BY DIVIDING الا‎ ۱۳۳۲ ۶ ٦ 
FUNCTION FRO ٠81ط‎ +4) ٣ 


DO 91 122 

DO 91 J=1,N2 
FMCF(1,J)=FMCF(1,J)/FNORM(I,J) 

CONTINUE 


DO 87 I=1,1 

DO 87 3-1 0 
en 

CONTINUE 


S 
(e 





FRENCH) 


( 6(6( 2 رہ 


رہ رہ رہ وہ رہ دہ 


DAD AD 


PAUSE 


MIS INN NS ES MUPTTHENÄRKÄYS TO PLOT THE 2-D MCF 


NTOT-2*NSIZE 
DO T11 
Bence, ٣۶ 
DS) 
FMAG(J)=FMCF(1,J) 
827 CONTINUE 


DO 83 J=1,NTOT 
WRITE(*,*)FMAG( J) ,DIST(J) 
83 CONTINUE 
PAUSE 


THE FOLLOWING SUBROUTINES ARE FOR THE GRAPHICS PACKAGE 
INOEWATMOSPHERTGTHEF IS BEING PLOTTED AT THIS POINT 


CALL VSINIT(18,8. ,10. ,0, ' MCF1. PLT' , IUNITV, IVID,5) 

CALL ORIGIN(.5,1. 5,0) 

CALL SCALE(FMAG,5. ,NTOT, 1) 

CALL AXIS(0.,0.,'MCF',0,1,1,5. ,90. ,FMAG(NTOT41), FHAG(NTOT*2) ,. 1,1) 
CALL SCALE(DIST,5. ,NTOT, 1) 


Pome cco 0.” ' 6.-1,-1,5. ,0. ,DISTCNTOT4 1), DISTCNTOT42) ,. 1,1) 
CALL LINES(DIST,FMAG,NTOT, 1, -1, ISYMB,. 1) 
CLOSE(IUNITV) 


CALL INT86(ITNO,IREG) 

CALL MSG(0.,0.,. 15,'PRESS ANY KEY TO CONTINUE',0.,0,1) 
ANS=GETC( ) 

CALL GMODE(IVID) 


THE MAGNITUDE OF THE FFT'D ELECTRIC FIELD IS CALCULATED IN ORDER 
MUNELOT THE OUTPUT. 


FMAX=0. 0 

DO 80 I=1,NR 

DO 80 J=1,NR 
FIELDM(I,J)=SORT(FIELDR(I,J)***2+FIELDI(I,J)**2) 
X-FIELDM(1,J) 
IF(X.GT. FMAX) THEN 
FMAX=X 
ENDIF 

80 CONTINUE 


THIS SECTION BEGINS THE CALLING SEOUENCE FOR PLOTTING 


CALL VSINIT(MON, 10. ,8. ,0, ' DITHER. PLT' , IUNITV, IVID,5) 
DO 100 I=1,N2 
DO 100 J=1,N2 

IX=J*INCR 

IY=I*INCR 

XMAX=ALOG10(FIELDM(I,J)/FMAX) 


N 
po 


С СО СО СО 


OOOO 


100 


275 


INDEX=8"ABS(XMAX)+1 

IF( INDEX. GE. 21) THEN 

ICOLOR=0 

ELSE 

ICOLOR=NDEX( INDEX) 

ENDIF 
CALL PIXEL( 1X, 17.120806» 
CONTINUE 
CALL MSG(0:a1 120155 1 ENON] 
CLOSE( IUNITV) 
CALL INT86(1TNO,IREG) 
CALL MSG(0.,0.,.15,'PRESS ANY KEY TO CONTINUE’ ,O. ,0,0) 
ANS=GETC( ) 
CALL GMODE(IVID) 


THIS IF STATEMENT QUES THE PROGRAM TO START OVER AGATi gs 
FRAUNHOFER OR THE CONVOLUTION FORM OF FRESNEL ARE USED 


IFC(NR GT. MODES HEN 
605108999 
ENDIE 


THIS IF STATEMENT QUES THE PROGRAM TO FULLY COMPLETE PROPAGATION 
THROUGH ALL THE SLABS IN THE TRANSFER FUNCTION FORM OF FRESENL 


1) 6٤٣٦ 
GO TO 1000 

ENDIF 

GO 190599) 

END 


SUBROUTINE QUADI(NR.PISDELMSH DXODY OWVIESDEMS 

COMMON /BLK2/ FIELDR( 256,256) ,ETELDE 258 ZEND 

DX-DELMSH 

DY-DELASH 

MID=(NR/2)+1 

00 27 +٤۶7 
2) ۲ - DX 

DOGS GI NE 
AI MID 8 8 DP 
THETA=PI*(C (X= X) ACY) (AVE DELLA) 
MX=FIELDRCU1 J) COS ITE 
YY=FIELDR( I J +٦٣ 
LA=FIELDICI II UOC ITE, 
WW=F IELDIC1I,J)*SIN( THETA) 
FIELDR(1,J)=XX-WW 
FIELDIN J)=YYTZZ 

CONTINUE 

RETURN 

END 


SUBROUTINE QUAD2(NR,PI,DX,DY,WVL,DELX) 

COMMON /BLK2/ FIELDR(256,256),FIELDI(256,256) 
DX22DX"WVL*DELX 

DY2=DY*WVL*DELX 

MID=NR/2+1 


tA 
[J 





CI CCI 


DO 274 I=1,NR 
Y=(I-MID)*DY2 

DO 274 J=1,NR 
X=(J-MID)*DX2 
PHI=PI*(((X*X)+(Y*Y)) /(WVL*DELX)) 
CX=FIELDR(1,J)*COS(PHI) 
CY-FIELDR(I,J)*SIN(PHI) 
CZ-FIELDI(I,J)*COS(PHI) 
CW-FIELDI(I,J)*SIN(PHI) 
CBR=CX-CW 
CBI=CY+0Z 
FIELDR(1,J)=(CBI/(WVL*DELX)) 
FIELDICI,J)=-1. *(CBR/(WVL*DELX)) 


274 CONTINUE 


RETURN 
END 


SUBROUTINE TRNSFR(NR,PI,DX,DY,WVL,DELX,DTRNS) 
COMMON /BLK2/ FIELDR(256,256),FIELDI(256,256) 
DX=DTRNS/NR 
DY=DTRNS/NR 
MID=NR/2+1 
DO 275 I-1,NR 
7۷7<) 1 7ھ‎ 
DO 275 J=1,NR 
FX=(J-MID)*DX 
FEE=-1. *PI"WVL*DELX*((FX*EX)+(EY*FY)) 
GX=FIELDR(I,J)*COS(FEE) 
GY-FIELDR(I,J)*SIN(FEE) 
GZ-FIELDI(I,J)*COS(FEE) 
GW=FIELDI(I,J)*SIN(FEE) 
FIELDR(1,J)-GX-GW 
FIELDI(I,J)=GY+GZ 
CONTINUE 
RETURN 
END 


SUBROUTINE FLTR(NR,DELMSH,CN2,DELX,WVL,FILTER) 
NOTE: DELMSH*NR IS THE LARGEST APERTURE SIZE 
COMMON /BLK3/ PHASER(256,256),PHASEI(256,256) 
PI=3. 141592653589792 

POWER=-11. /6. 
TPI=2. *PI 

N2=NR/2 

NPIVOT=N2+1 

LAST=NPIVOT+1 

DLKAPA=(TPI/(NR*DELMSH) )**POWER 
FACTOR=SORT((TPI***:3)*, 033**CN2*DELX/(WVL**2)) 
FUDGE=DLKAPA*FACTOR 

DO 100 I=1,NPIVOT 

EYE=REAL(I) 

EYE2-EYE*EYE 

DO 100 J=1,NR 

IF(J.LE.NPIVOT) THEN 
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WHY=REAL(J) 
ELSE 
WHY=REAL(NR-J+2) 
ENDIF 
XKAPPA=SQRT(EYE 2+WHY*WHY) 
AKAPPA=XKAPPA"*FILTER 
PHASER(1,J)=PHASER( 1,J)*FUDGE*AKAPPA 
PHASEI(1,J)=PHASEI(1,J)*FUDGE*AKAPPA 
100 CONTINUE 
DO 110 I=LAST, NR 
EYE=REAL(NR-I+2) 
EYE2-EYE*EYE 
DO 110 J=1,NR 
IF(J. LE. NPIVOT) THEN 
WHY=REAL(J) 
ELSE 
WHY=REAL(NR-J+2) 
ENDIF 
XKAPPA=SQRT( EYE2+WHY*WHY ) 
AKAPPA=XKAPPA***F ILTER 
PHASER(1,J)=PHASER(1,J)**FUDGE*AKAPPA 
PHASEI(1,J)=PHASEI(1I,J)*FUDGE**AKAPPA 
110 CONTINUE 


PHASER(1,1)=0.0 
PHASEI(1,1)=0.0 


NOTE: .03441=TPI**-11./6. TO TRANSFORM FROM KAPPA TO FREQ SPACE. 
DO 11 I=1,NR 
DO 11 J=1,NR 
PHASER(1,J)=PHASER(1,J)*. 03441 
PHASEI(1,J)=PHASEI(1,J)*. 03441 
11 CONTINUE 
RETURN 
END 


SUBROUTINE IFTSCR(NR,M,DELMSH) 
COMMON /BLK1/ RE(256),RIM(256) 
COMMON /BLK3/ PHASER(256,256),PHASEI(256,256) 
PI=3. 141592653589792 
TPI=2. *PI 
Sie) eno 
DO 20 I=1,NR 
DO 21 J=1,NR 
RE(J)-PHASER(1,J) 
RIM(J)=PHASEI(1,J) 
21 CONTINUE 
CALL FFT(M,SIGN,DELMSH) 
DO 22 J=1,NR 
PHASER(1,J)=RE(J) 
PHASEI(I,J)=RIM(J) 
22 CONTINUE 
20 CONTINUE 
DO 30 J=1,NR 
DO 31 T=1,NR 
RE(I)-PHASER(I,J) 


UA 
Eja 





DT 


352 
30 


51 


62 
60 


71 


2 
70 


RIM(I)=PHASEI(1,J) 
CONTINUE 
CALL FFT(M,SIGN,DELMSH) 
DO 32 I=1,NR 
PHASER(1,J)=RE(1) 
PHASEI(I,J)=RIM(I) 
CONTINUE 
CONTINUE 
RETURN 
END 


SUBROUTINE DFTIFT(NR,M,SIGN,DELMSH) 

COMMON /BLK1/ RE(256),RIM(256) 

COMMON /BLK2/ FIELDR(256,256),FIELDI(256,256) 

DO 60 I=1,NR 
DO 61 J=1,NR 

RE(J)-FIELDR(1,J) 
RIM(J)=FIELDI(1,J) 

CONTINUE 

CALL FFT(M,SIGN,DELMSH) 

DO 62 J=1,NR 
FIELDR(I,J)=RE(J) 
FIELDI(I,J)=RIM(J) 

CONTINUE 

CONTINUE 

DO 70 J=1,NR 
DO 71 I=1,NR 

RE(I)-FIELDR(1,J) 
RIM(I)-FIELDI(1,J) 

CONTINUE 

CALL FFT(M,SIGN,DELMSH) 

DO 72 I-1,NR 
FIELDR(I,J)=RE(I) 
FIELDI(I,J)-RIM(I) 

CONTINUE 

CONTINUE 

RETURN 

END 


SMEROUTINE FFT(M,SIGN,DELMSH) 
BOJNOM /BLKI/ RE(Ż56),RIM(256) 
eo. 141592653589792*SIGN 


0 71 
шее 
3-1 


DO 200 I=1,N1 
agent. J) THEN 

T=RE(J) 
RE(J)=RE(I) 
RE(1)=T 
T=RIM(J) 
RIM(J)=RIM(I) 
RIM(I)=T 

END IF 

K=N/2 


201 


200 


204 


207 
202 


ho 
O 
Ui! 


301 


300 


DO 201 WHILECK. LT. J) 
J=J-K 
PH 
CONTINUE 
J=J+K 
CONTINUE 
LES! 
DO OT 
LEl=TE 
LE=LE+LE 
URE-1. 
UIM=0. 
A CEPI 
WRE=COS(ANG) 
WIM=SIN(ANG) 
DO 203 J=1,LE1 
DO 204 I=J,N,LE 
IP=I+LE1 
TRE-RE(IP)*URE-RIM(IP)*UIM 
TIM=RE(IP)*UIM+RIM(IP)*URE 
RE(IP)-RE(I)-TRE 
RIM(IP)=RIM(I)-TIM 
RE(I)=RE(I)+TRE 
RIM(I)=RIM(I)+TIM 
CONTINUE 
T=URE*WRE -UIM*WIM 
UIM=URE*WIM+UIM*WRE 
URE=T 
CONTINUE 
CONTINUE 
1۳۷ ۰۱۳ ۰ 
PTS=1. 0/(N**DELMSH) 
DO 205 I=1,N 
RECI)-RECI)*PTS 
RIM(I)=RIM(I)*PTS 
CONTINUE 
ENDIF 
RETURN 
END 


SUBROUTINE GGAUS(NR,SEED) 
COMMON /BLK3/ PHASER(256,256),PHASEI(256,256) 
DO 300 I=1,NR 

DO 300 J=1,NR 

V1=2. *RAN(SEED)-1 

V2=2. *RAN(SEED)-1 
S=V1*V1+V2:V2 

۷ ٰ9۹پ٤‎ corner 201] 
SCALE=SQRT( -2. *ALOG(S)/S) 
X1=/1*SCALE 
X2-V2*SCALE 
PHASER(1,J)=X1 
PHASEI(1I,J)=0. 0 

CONTINUE 

RETURN 
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END 


SUBROUTINE MCF(FNORM,NR,M,PI,NSIZE,N2,DELMSH) 
CONNOR UBER БЕ( 2556 К КТ1М( 256) 

COMMON БЕКУ БЕ ПЕШКЕ 2567250), FIELDI( 256,256) 
DIMENSION "ENORM( 256) 256) 


MISS EEIIONCREAPESNLHE PLANAR ELECTRIC FIELD 


DO 39 I=1,NR 
DO 39 J=1,NR 
FIELDR(1,J)20.0 
FIELDI(I,J)=0.0 
CONTINUE 


DO 45 I=N2-NSIZE+1,N2+NSIZE 

DO 45 J=N2-NSIZE+1,N2+NSIZE 
ETELDRO J)=130 

CONTINUE 


SIGN=-1.0 
CALL DFTIFT(NR,M,SIGN,DELMSH) 


DO 80 I=1,NR 

DO 80 J=1,NR 

ENORMCI ,J)=FIELDRCI,J)**2+FIELDI(I,J)**? 
CONTINUE 


DO 123 I=1,NR 
DO 123 J=1,NR 

FIELDI(1,J)20.0 
CONTINUE 


SIGN=+1. 0 

DO 90 I=1,NR 

DO 91 J=1,NR 
RE(J)EFNORM(1,J) 
RIM(J)=FIELDI(1,J) 

CONTINUE 

CALL FFT(M,SIGN,DELMSH) 

DO 92 J=1,NR 
FNORM(1,J)=RE(J) 
FIELDI(1,J)=RIM(J) 


CONTINUE 

CONTINUE 

DO 93 J=1,NR 

DO 94 I=1,NR 
RE(1)=FNORM(I,J) 
RIM(I)-FIELDI(1,J) 

CONTINUE 


CALL FFT(M,SIGN,DELMSH) 
DO 95 I=1,NR 
FNORM(I, J)-RE(I) 
FIELDI(1,J)SRIM(I) 


Uh 
~J 
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CONTINUE 
CONTINUE 


FLDY=0. 0 
DOR NR 

DO 88 J=1,NR 

XMG=FNORM(I,J) 

IF(XMG.GT. FLDM) THEN 

FLDM=XMG 
ENDIF 

CONTINUE 
DO 89 I=1,N2 

DO 89 J=1,N2 
FNORM(I,J)sFNORM(I, J)/FLDM 
CONTINUE 


RETURN 
END 
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